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0.  ABSTRACT 

i  ff  ^ 

Multidimensional  Wilson- Hilfertyq32f  transformations  support  Gaussian  approxi¬ 
mations  to  certain  joint  distributions  of  quadratic  forms  in  Gaussian  variables.  Central  and 
noncentral  distributions  are  studied  and  applications  are  noted.  Parameters  of  the  approx¬ 
imating  distributions  are  given  up  to  terms  of  order  f(v~3)  in  the  degrees  of  freedom  v.  v 
Numerical  studies  validate  using  these  approximation^  over  a  range  of  parameters  for  an 

essential  subclass  of  the  distributions.  ^  _ , 

Orv,'*-.i~v 

I.  INTRODUCTION 

Developments  in  statistics  and  applied  probability  often  entail  joint  distributions  of 
definite  quadratic  forms  in  Gaussian  variables,  either  m  small  samples  or  asymptotically. 
Examples  include  linear  statistical  models  ( 10].  the  ballistics  of  multiple  weapons  systems 
|4],  signal  detection  in  multichannel  receivers  [25],  and  bone  iengths  determined  in  vivo  using 
X-ray  stereography  ((26|,  (2SJ  and  (29|).  Topics  in  large-sample  theory  include  limiting  joint 
distributions  of  likelihood  ratios,  of  Pearson  s  [27]  A7  statistics  for  categorical  data  |12],  and 
of  Friedman  s  (5]  statistics  in  two-way  data  without  normality  (|  13]  and  1 14]).  Further 
details  are  given  subsequently,  and  other  examples  could  be  cited. 

Distributions  of  these  types  include  multidimensional  chi-squared  and  Rayleigh  dis¬ 
tributions  ([8],  |9],  [10],  [18],  [19],  [20],  [22|,  [24],  and  [30]),  often  depending  on  an  excess 
of  parameters.  Series  expansions  for  their  distribution  functions  typically  are  intractable  in 
dimensions  greater  than  two:  convergence  properties  of  these  series  may  not  be  known  (19); 
and  series  known  to  converge  may  do  so  slowly.  Viable  approximations  to  these  distrib¬ 
utions  are  clearly  needed.  Among  alternatives,  multivariate  Edgeworth  expansions  [3]  are 
often  flawed  by  inadequacy  of  the  leading  term  and  failure  of  the  sum  of  the  fust  few  terms 
to  be  positive.  Approximations  from  the  Pearson  system,  in  wide  use  in  the  univariate 
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case,  hold  scant  promise  since  little  is  known  beyond  the  bivariate  case  ([18],  pp.  6-9). 
Alternatively,  the  normalizing  transformations  of  Wilson  and  Hilferty  [32]  appear  promis¬ 
ing,  as  these  give  remarkable  accuracy  for  a  single  quadratic  form  over  a  wide  range  of  pa¬ 
rameters  ([16],  [17]  and  [23]). 

Here  we  develop  multivariate  transformations  leading  to  Gaussian  approximations 
not  depending  on  excessive  parameters,  requiring  only  moments  of  first  and  second  orders. 
To  be  precise,  let  [£?,, ... ,  Qr]  be  standardized  quadratic  forms  in  Gaussian  variables  having 
v  degrees  of  freedom  whose  limit  as  v  -+  oo  is  jointly  Gaussian.  Lemma  1  assures  that 
IQI\  ... ,  <2’r]  also  has  a  Gaussian  limit  for  every  {a,  e  (0,1];  1  <j  <  r).  Extending  devel¬ 
opments  in  [16],  we  consider  multivariate  power  transformations  {£?,  -» (<2,)“J;  1  <  y  £1  r} , 
choosing  («,, ... ,  a,}  so  as  to  accelerate  the  convergence  of  certain  moment  sequences  to 
those  of  Gaussian  variables.  Parameters  of  the  approximating  distributions  are  given  up  to 
terms  of  order  o(v3),  and  the  approximation  is  studied  numerically  in  selected  cases. 


2.  PRELIMINARIES 


2.1  Notation.  Standard  spaces  include  the  Euclidean  n-space  R",  the  space  F„,p  of  real 
(n  x  p)  matrices,  the  symmetric  (p  x  p)  matrices  Sp,  and  the  cone  Sp  <=■  Sp  of  positive  semi- 
definite  varieties.  The  half-open  unit  cube  in  Rn  is  denoted  by  D(n)  = 
(0,  1]  x  (0,  1]  x  ...  x  (0,  1].  Matrix  operations  include  transposition  ( A  -*  A')  and  inversion 
{B-*  B-').  Special  arrays  are  /„,  0,  A  x  B,  and  Diag(T,, ...  ,Ar),  denoting  respectively  the 
( n  x  n)  identity,  a  matrix  of  zeros,  the  direct  product  [a,, 5],  and  a  block-diagonal  matrix, 
in  addition  to  the  unit  vector  l„  =  [l, ... ,  1]'  e  R*. 

2.2  Some  Basic  Distributions.  JZ(Z)  denotes  the  law  of  distribution  of  Z.  The  joint 

moments  of  X  e  R'’  are  denoted  by  ,  its  central  moments  by  and  its  joint 

cumulants  bv  k;i  , ,  all  of  order  s  =  s,  -i - (-  s„.  We  consider  the  probability  density  func- 

tion  (pdf),  the  cumulative  distribution  function  (cdf),  and  the  cumulant  generating  function 
(cgf)  of  Z.  Gaussian  laws  on  R'  and  Fn ,,  are  denoted  by  Np(p,  £)  and  A 'n,,(.V/,  D>  respec¬ 
tively,  where  p  e  R'  and  M  e  F„,p  are  arrays  of  means,  and  £  e  Sp  and  T  e  5^,  are  arrays  of 
dispersion  parameters.  The  Wishart  distribution  on  Sp,  having  v  degrees  of  freedom  and 
the  matrices  £  and  ®  of  scale  and  noncentrality  parameters  in  Sp,  is  denoted  by 
H^(v,  £,  0).  In  particular,  if  St(Y)  =  Ar„,(/W, L  x  I),  then  <£ (E'E)  =  W^n,  ?.  WM)- 


For 


2.3  Quadratic  Forms.  We  represent  quadratic  forms  in  Gaussian  variables  via  parti-  itl 
tioned  Wishart  matrices.  Suppose  that  y(y)  =  Np(p,Z);  consider  quadratic  forms 
\y'A$,...,y'A,y\  such  that  {T,  e  Sp\  \  <j<  r}\  identify  the  joint  distribution  of 
z-  [r,\ ...  ,2/7  on  P/p  as  Nrp(0,  T),  where  {zt  =  AJl2y;  \  <j<r),  0  =  [0, 0/]'  and 
{0,  =  AJ,2p;l  <j<r}\  and  write  T  e  5 *  in  partitioned  form  as  T  =  [.4,I/J£.T'3].  If 


od 
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2  llty  Codes 
.1  and/or 
Special 


W  =  zf  %  then  ^ity)  =  Wrp(l,r,00').  Partitioning  H7  =  [  H7’,]  conformably  with 
[2,', ... ,  2/]'  shows  that  [y'Aj, ...  ,y'Ary]  may  be  represented  equivalently  as 
LSiZi . z/zj  and  [trip,, . trJT,,]. 

Generally  [trJTn, ... ,  trW7,,]  are  quadratic  forms  in  jointly  Gaussian  variables.  For  if 
W  =  Y'Y  with  S£[Y)  =  A L  x  Z),  then  ^(H7)  =  W7,(/i,  I,  ©)  with  ©  =  M'M.  Par¬ 
titioning  the  typical  row  of  Y  as  y\  =  [y',  ...  ,y',  ]  such  that  ytj  e  R'’-'  with  p,  +  -  +  pr=  p, 

and  thus  Y  =  [K, ...  >',]  with  {T,  e  Fn, . ;  1  <;  <  r},  we  infer  that  { IIP  =  ;  1  <j  <  r), 

*  J  i=r  ~ 

so  that  {trlf7;  =  1  < y  S  r)  are  quadratic  forms  in  jointly  Gaussian  variables. 

i=r 

More  general  forms  {tr!j77Z?y;  1  <j  <  r),  with  (5,  e  S;;  1  <  y  <:  r} ,  are  also  quadratic 
forms  in  jointly  Gaussian  variables  because  {trW7,#,  =  5, >>,■,;  1  <  y  <  r} .  The  trans¬ 

formation  Y-*YB',  with  B'  =  Diag(5j'2, ... ,  B}12),  shows  that  YB'  has  the  distribution 
S£(Yff)  =  Nn,,(MB\  i,  x  BLB'),  and  thus  K=Z'Z,  with  Z=YB',  has  <£{V)  = 
Wf(n,  BLB',  B@B').  Partitioning  V—[VvJ  conformably  with  H'=  [W7y]  shows  that 
{trV,  =  trWjjBj',  1  <  j  <,  r}  are  positive  semidefmite  quadratic  forms  in  Gaussian  variables, 
so  that  {/?!, ... ,  B,j  may  be  absorbed  into  parameters  of  the  Wishart  distribution.  In  partic¬ 
ular,  the  case  n  =  1  yields  the  forms  {}\'B]yu  ...  ,yr'B,yr~\  iny  =  [y\ . y\  ]'  on  R"  as  before. 

Representations  in  all  these  cases  can  be  given  in  terms  of  standardized  variables. 
Specifically,  the  quadratic  forms  [y'A^y, ...  ,y'A,y ]  in  y  e  R'\  with  if(y)  =  Nf(n,  Z),  become 
[z’SiZ, ... ,  ?B,z]  with  {. B ,  =  Vl2A£'i2;  1  <j  <,  r}  and  g  =  I-1'2;*  such  that  i?(z)  =  A’,(0,  /,) 
with  6  —  The  further  transformation  z  -+  P/z  ~  where  Pt  is  an  orthogonal  matrix 

chosen  to  diagonalize  Z XI2AJL'12  -*  PI%'[I2AJLXI2P]  =  A,  =  Diag(dyl  ...  ,  6JP),  yields  equivalent 
forms 

{/Ajy  ~  Z'BjZ  =  u'j  A jU/,  1  <j<r)  (2.1) 

in  the  Gaussian  vector  u  =  [u7,  , ... ,  t/r  ]'  on  Rr?.  Here  i?(u)  =  A'rf(y,  Z)  with  y  = 

[y; . y',1'  such  that  { y,  =  P;'Z_,,2fx;  1  <y  ^  r),  and  in  partitioned  form  E  = 

[P/jP,]  e  Other  cases  may  be  treated  similarly. 

2.4  A  Limit  Rule.  The  following  lemma  is  basic,  where  <1 >,(z)  denotes  the 
p-dimensional  Gaussian  cdf  with  moments  to  be  specified,  and  Z'  denotes  the  random  di¬ 
agonal  matrix  Z’  =  Diag  {Z\\  ... ,  Z‘/)- 

LEMMA  L  Let  {ZN\  A' =  1,2,...}  be  a  sequence  of  random  diagonal  (pxp)  matrices 
having  nonzero  means  E{Zh)  =  MN  and  the  diagonal  matrix  Z;v  containing  variances  only, 
and  let  PA.  =  z.~J{ZN-  MN)  be  such  that  (i)  lim  =  C)(y)  and  (ii)  lim  A/a’Ev=  0.  If 

=  Diag(Z“’, ... ,  Z*'),  then  for  eve  ’  {a.  e  (0, 11;  l^y  <p}  we  have  lim  if(lf'A<2)) 

¥  *  S —00  ~~ 

=  ^(b7). 
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Proof,  Let  A  =  Diag  (a, . a,);  write  Zv  =  EvFy  +  A/v;  and  note  on  dropping  subscripts 

that  Z*  =  A/*(7  +  A/-,Z)')'.  Expanding  each  element  of  (/  +  M-'E.Y)*  in  its  binomial  series, 
using  hypothesis  (ii)  for  N  sufficiently  large  and  the  assumption  that  a.J  e  (0,  1],  gives 

<7  +  A/_1Sy)“  = /  +  AM~%Y  +  o(\r]-EY).  (2.2) 

On  combining  expressions  and  reinstating  subscripts,  we  have 

iXsr  —  A/y)  —  dX!s  HnXn  "b  N  (2-3) 

The  lemma  now  follows  from  hypothesis  (i)  on  taking  limits  in  (2.3),  where  it  is  seen  that 
Z;  asymptotically  is  linear  in  the  limiting  Gaussian  matrix  and  where  an  additional 
standardization  may  be  used  if  needed  to  secure  a  proper  limit  for 


3.  NORMALIZING  TRANSFORMATIONS 


We  seek  normalizing  transformations  [Q;  — *  1  <j  <r)  for  standardized  quad¬ 

ratic  forms  in  jointly  Gaussian  variables.  From  Section  2.3  it  suffices  to  consider 
[t/,/v0|, ... ,  UJvd,'),  where  {L,  =  trlj'y  1  <j<r],  such  that  f£(W)  =  W^v,  I,  ©)  and 
{£(trW7,)  =  v0.;  1  <  y  <  r] .  We  first  evaluate  joint  cumulants  of  [  (/, , ... ,  LJ  from  their  cgf; 
we  next  convert  cumulants  to  moments;  and  we  then  expand  typical  moments  of 
{Vj  =  (i/,/v0;)*';  1  <j  <  r}  to  terms  of  required  order  in  v.  We  frnally  choose  {a,, ... ,  a,}  so 
as  to  accelerate  the  convergence  of  certain  moment  sequences  to  those  of  jointly  Gaussian 
variables.  We  outline  these  steps  in  the  sections  following,  deferring  further  details  to  Ap¬ 
pendix  A. 

3.1  The  Joint  Moments.  Starting  with  the  joint  chf  of  W=  [H7,]  with  argument 
T  e  Sp  ( cf.  |2])  and  expanding  its  cgf  as  in  1 1 1).  we  write  the  joint  cgf  of  [(/,, ... ,  with 
argument  t  =  [/,, ... ,  /,]',  as 

^ £/({)  =  v  I  /Vr'trf  (7'05)J-'j0(Z  +  50)]  (3.1) 

s=\ 


where  Q  =  v_)©  and  T0  =  Diag  (tx[p . .  lr/p).  The  joint  cumulant  x  ,r,  of  order  s  = 

r,  H f-  sr,  is  the  coefficient  of  i’f' ...  tfjsf. ...  r,!  in  the  power  series  (3.1).  To  simplify  no¬ 
tation,  partition  W  =  [W7,],  I  =  [I,;]  and  0  =  [0,;]  conformably;  let  9,  =  tr(2„  +  Q;;)  and 
o)jj  =  tr +  sQjj)  for  1  <  y  <  r,  define 


°s(h‘2 hh  -  v'i)  =  Xrhi2*hh  -  +  ^y,) 


(3.2) 
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where  s' =  s  —  1  and  s=  1,2,...;  and  write  o^ijjiit)  =  a^ifp),  for  example.  The  joint 
cumulants  kS[  ...  v  of  [L’,t ...  ,  of  orders  1  <  s  <  4,  are  summarized  in  Table  1  using  the 
symmetry  of  Z  and  Q  as  in  Appendix  A. 

TABLE  1.  Typical  cumulants  kS|  ,r  of  [t/, . L'J  of  order  j,  1  <  s  <  4. 

r  /Ji ...  Hr  Coefficient  of  i‘t\\  ...  r^/j,! ...  jr! 


1 

r, 

v6i 

2 

1} 

2va >a 

2 

2v<rj  (ry2/) 

3 

8vO>(3 

3 

8v[«3(/y»P)  +  2<7J(i3/Ji)]/3 

3 

8v[o3(//TW)  +  <T3(i*2/20  +  ^3W)]/3 

4 

if 

48vo>,4 

4 

24v[<r4(/5y2/)  +  ct(PfP)2 

4 

w 

8  v[2<t4(£//)  +  2o4(ifPj*i)  +  <r4(ij*P)  +  <r4(/?/)] 

4 

8v[o4(/^2A  +  ai(Pk1pI)  +  at(ij2Pk2i) 

-(-  at(ij7k2P)  4-  oa(jPk2j )  4-  oa(ik2Pj2i)} 

Moments  of  [O',, ... ,  0,]  are  computed  from  Table  1  using  known  relations  between 
moments  and  cumulants  as  in  Appendix  Table  Al;  see  [6],  for  example.  These  computa¬ 
tions  give  typical  joint  moments  of  orders  1  <  s  <,  4  as  summarized  in  Table  2,  and  other 
moments  follow  similarly. 


TABLE  2.  Typical  joint  moments  of  [(/,,  U2,  C/}]  of  order  r,  1  <  s  <  4. 


3.2  The  Transformations.  In  what  follows  suppose  that  0  =  o(v)  as  v-*oo,  and 
conside  he  standardized  variables  [L',/v0, . 6',/v0,]  with  {v0;  =  E{Vj)\  1  <j  <  r}  to  en- 


sure  a  proper  limit.  We  establish  in  Theorem  1  that  [(L^/v#,)’1, ...  ,(Urlvdr)','\  are 
asymptotically  Gaussian  for  every  a  e  D{r). 

THEOREM  !.  Suppose  that  !£{W)  =  B'(v,  L,  0)  such  that  0  =  o(y )  as  v  -+  oo.  Then  the 
variables  [Lr,/v0,,...  ,  UJvd,'],  with  {Uj  =  \xWn\  1  <j<r},  are  asymptotically  Gaussian  as 
v  oo,  as  are  [(t/j/vfl,)*1, ...  ,(t/,/v0,)*']  for  every  a  e  D(r). 

Proof.  The  variables  W  =  Y'Y  when  standardized,  and  thus  [ UJvd „ ... ,  UJv9r]  as  mar¬ 
ginals,  are  asymptotically  Gaussian  by  central  limit  theory  [11],  verifying  condition  (i)  of 
Lemma  1  for  Z,  =  Diag  (Uxlv9u  ... ,  L',/v0,).  To  verify  condition  (ii),  identify  A/,  and  Zv 
as  in  Lemma  1  with  N  =  v  and  observe  from  Table  1  that 


hm  =vhm  (2/v)’/2Diag  (y, . y,)  =  0  (3.3) 

where  {y,  =  a>)J2  Bjx\  1  <j  <r)  and  y  =  o(l)  because  0  =  o(v).  Lemma  1  now  asserts  that 
the  transformations  £(«):{  L)/v0,  -» (bjvdy>\  1  <j<r)  yield  asymptotically  Gaussian  vari¬ 
ables  for  every  {«;  e  (0,1];  1  <j  <  r),  which  completes  our  proof. 


It  remains  to  examine  certain  moment  sequences  of  [( L/,/v0,)*1 . (i/,/v0,)*'],  which 

we  now  investigate.  Letting  { bt  =  v0y  -f  ey;  \  <j  <.r]  and  { V/  =  (t/y/vfl;)*>;  \  <j<r j,  we 
evaluate  the  mixed  moments 


E{v\l ...  Vs;)  =  £  [(1  +  e./ve,)”’1’  ...  (1  +  erl*er)*' Srl  (3.4) 

on  expanding  each  factor  and  finding  expected  values  term-wise.  For  moments  to  third 
order,  it  suffices  to  consider  the  case  r  =  3  and  to  use  {a,  /?,  y)  in  lieu  of  {a,,  a  2,  a}).  The 
resulting  expression  is  basic,  namely, 


E(K^3)  =  I  S  I  ( T ^eVeJe^albjclEie^) 

1=0  j=Q  k=0 


(3.5) 


where  {a',  bj,  cl)  are  binomial  coefficients  identified  in  Appendix  A. 
TABLE  3.  Marginal  moments  of  Vx  =  (trl Yujv9xy  to  terms  of  order  o(v-3). 


Moment 

Expression  * 

Mioo  («) 

l+  ,8,  (“‘"'W  +  .’8j  [(  3  )*’+(  2_)«] 

Mjoo(«) 

[2*,  +  3(«-l> W 

**[,]  =  «(«  -  i)  -  (« -  0; 

4>i  —  '(?n  +  jQn)/tr(Ln  +  On). 
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To  continue,  identity  moments  about  zero  as  n'„( a, /?,>•)  =  let 

Mr„(a,  /?,  y)  be  the  corresponding  central  moment;  and  let  nrs0 (a,  /?),  /*„,,(/?,  y)  and  m*>o(») 
typical  marginal  moments.  We  next  compute  (a,  /?,  y);  r  +  s  +  t  =  1}  and 
{/*„,( a,  /?,  y);  2  <  r  4-  j  +  t  <  3}  to  terms  of  order  o(v-3),  applying  expression  (3.5)  as  often 
as  needed  to  known  relations  between  central  and  noncentral  moments  and  using  entries 
from  Table  2.  Moments  of  each  one-dimensional  marginal  distribution  are  taken  from  ( 1 6] 
on  replacing  their  0,  by  our  v0;,  for  example.  With  at0  =  a(a  —  1) ...  (a  —  l)  and 
{4>,  =  (o,Jv6 s  =  2,  3, ... },  these  moments  are  given  in  Table  3  to  terms  of  order  o(v3)  for 
(t/,/v0,)*,  with  corresponding  expressions  for  (t/2/v02/  and  (f/3/v03)?. 


TABLE  4.  Nonvanishing  terms  to  order  i  +  j<L  4  and  their  coefficients  in  the  expression 
(A  .6)  for  the  joint  moment  ^i,10(a,  /?)  of  [L,,  V2,  Lj], 


‘  j 

<+>■ 

+2  (0;0i)-‘ 

a)b } 

r — i 

i 

i — i 

1  1 

v-J 

(W 

*P 

2v<t2(  1 22 1) 

2  1 

V-3 

(«-• 

«(<*  -  l)/?/2 

8v[<t3(12213)  +  2a3(l322l)]/3 

1  2 

V-3 

(0,01)-' 

*P(P  ~  l)/2 

8v[cr3(2 1223)  +  2ct3(23122)]/3 

2  2 

v~* 

(02,oi)" 

«(«  -  W(P  ~  D/4 

8v[^:o  +  va22(122l)] 

3  1 

V-4 

(d'A)-' 

oc(a  -  1)(«  -  2)0/6 

12v[2M310  +  vw12ct2(1221)] 

1  3 

V-4 

(0,01)" 

8 

/**s 

*C& 

1 

1 

ON 

12v[2A'130  4-  vcoJ2<t2(1221)] 

TABLE  5.  Nonvanishing  terms  to  order  i+j <>  4  and  their  coefficients  in  the  expression 
(A.l)  for  the  joint  moment  Mjio(a>  P)  °f  f'j.  L3]. 


h 

i 

j 

(_L 

(8T0$" 

Coefficient 

[£(ef)£(e ;ei)  -  £(ef)£(e;)£(e£)3 

0 

1 

1 

v-2 

(0A)" 

0 

2va2(122l) 

0 

2 

1 

V-3 

(02A)" 

a2/? 

Sv[tr3(122l3)  4-  2a3(  l322l)]/3 

0 

1 

2 

V-3 

(0,01)" 

0 

8v[u3(21223)  4-  2ct3(23122)]/3 

0 

2 

2 

v-4 

(0]0D" 

-  l)/2 

8v[A220  +  v<t2(1221)] 

0 

3 

1 

V-4 

(0\0:)" 

0a2(a  —  1) 

12v[2A'3)0  +  va>,2a2(122l)] 

0 

1 

3 

V-4 

(0,0\)" 

0 

12v[2A)30  +  va>22er2(122l)] 

2 

1 

3 

V-4 

(0\O2)-' 

1 

8 

K» 

1 

4v2co12a2(122l) 

We  summarize  in  Tables  4.  5  and  6  the  computations  needed  for  Muo(a>  P)< 
H7l0(a,  p),  and  finl(a,p,y)  up  to  terms  of  order  o(v'3),  where  expressions  for  are 

identified  in  Table  7.  Details  are  given  in  Appendix  A.  The  required  partial  sums  for 
P,  y)}  are  found  from  each  table  on  ignoring  the  first  two  columns,  next  multiplying 
the  remaining  expressions  in  each  row,  then  summing  these  products  over  rows,  and  finally 
collecting  terms  in  powers  of  v1. 
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TABLE  6.  Nonvanishing  terms  to  order  i  +  j  +  k  <  A  and  their  coefficients  in  the  ex¬ 
pression  M.8)  for  the  joint  moment  /im(a,  P,  y)  of  [F,.  F2,  F3] 


i  j  k 

a)b)c\ 

C(i,j,  k)* 

1  1  1 

v-3  (Wj)'1 

aPy 

8v[a3(12232l)  +  <t3(  1 3222 1 )  +  <r3(212322)]/3 

2  1  1 

v-  (WlL1 

«(«  -  Wy/2 

8v[A'2II  +  vct2(1221)<t2(13j1)]** 

121 

7 

<—> 

7 

> 

ap(p  -  l)y/2 

8v[A)21  +  v<t2(21j2)<t2(23:2)] 

1  12 

v-  (0AGI)-' 

aPy(y  -  l)/2 

8v[  AT,|2  +  v<72(3123)a2(3223)] 

k)  is  defmed  in  expression  (A. 9). 

"This  expression  derives  from  {8vA2ll  +  4v2co12(t2(2322)  +  8v2<t2(  1 22  1)<t2(  1 32 1)  —  0  — 
(2vco12)[2v<72(2322)]  —  0  +  0} ,  for  example. 


TABLE  7.  Identification  of  typical  constants  {K,lk}. 


i 

j 

k 

3 

1 

0 

<y4(  ls22 1)  +  <t4(132213) 

2 

2 

0 

2cx4(13241)  +  2<t4(12212221)  +  a^l3)  4-  <r4(21«23) 

2 

1 

1 

[a4(l32232l)  +  <j„(l33222l)  +  a4(  122l232l) 

+  <r4(I2232l3)  +  a4(2  F322)  +  ct4(13212221)] 

3.3  The  Approximating  Distribution.  It  remains  to  choose  {a,/?,y}.  As  all  mixed 
central  moments  of  [F,,  F2,  F3]  of  order  3  vanish  up  to  terms  of  order  o(v~2)  (compare 
Tables  5  and  6),  we  specifically  choose  {a,  p,  y)  so  as  to  annihilate  the  leading  term  in  ex¬ 
pressions  for  each  of  nm(a.),  Mojo(^)>  and  nKi{y).  The  solution  from  Table  3  is 
a  =  1  —  201w13/3cu?2,  with  similar  expressions  for  P  and  y.  For  determining  {<%,, ...  ,  ar}  in 
the  general  case,  identical  arguments  yield 

=  1  —  26j(i)jjl3u>j2',  1  <j  <  r.  (3.6) 

The  foregoing  developments  support  a  Gaussian  approximation  to  the  joint  distrib¬ 
ution  of  [F,, ...  ,  F,].  Up  to  terms  of  order  o(vJ),  the  means  n  =  ... ,  /*,]'  and  dispersion 

parameters  Z  =  of  the  approximating  distribution  are  given  in  the  following  theorem; 
expressions  to  order  o(v-3)  are  available  from  Table  3  and  Table  4. 

THEOREM  2.  To  terms  of  order  o(vJ),  parameters  of  the  Gaussian  approximation  to  the 
joint  distribution  of  [F„ ... ,  F,],  with  {F;  =  (tr  WJvQfi-,  1  <j  <  r),  are  m  =  [ju„  ... ,  Mr]' 
and  Z  =  [£J,  where 


Hj  =  \  +  -  l)/v0?;  I  <,j<,r 

(3.7) 

Zij  =  2a/  aj  o2(i/20/v^;  i  <  h  j  r 

(3.8) 
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for  every  {*,  e  (0,  1];  1  <  j  <r). 

In  summary,  the  joint  distribution  of  [I7,, ... ,  Vr]  is  approximately  Gaussian  as  indi¬ 
cated.  Equivalently,  the  variables 

Zj  =  v]l26j[{UjlvQj)  1  -  1  -  a)j2a.j(aj  -  l)/v6y2]/(2a)y-2  a y)1/2;  1  <j  <  r  (3.9) 

may  be  taken  to  be  approximately  Gaussian  having  zero  means,  unit  variances,  and  the 
correlation  matrix  R  =  [pv],  where 

Pij  =  ^  2  ( (/' 2  0  /  co  2  coj^2 1  <i<j<r  (3.10) 

from  p/j  =  tiJV.PVJ7  and  (3.8).  To  terms  of  order  o(v~2),  these  correlations  do  not  depend 
on  the  particular  choice  for  {a,, ... ,  a,}. 

3.4  Some  Special  Cases.  Special  distributions  merit  further  attention.  As  in  Section 
2.3,  partition  W  =  [WA],  ?  =  [?,,],  and  ©  =  hi' M  =  [0,J  conformably;  let 

{}.j  =  tr©„;  1  <j  <  r};  and  suppose  that  (X;;  =  [p  ;  1  <j<r).  Then  the  marginal  distrib¬ 
ution  of  Uj  =  trWA  is  noncentral  chi-squared  having  vpj  degrees  of  freedom  and  the 
noncentrality  parameter  //,  the  corresponding  central  distribution  of  [trW7,,, ... ,  trlT),]  has 
been  studied  in  [10].  We  determine  from  (3.6)  that 

{<xj  =  1  -  2 {pj  +  +  UMlipj  +  2 ?.j)\  1  <y'<r).  (3.1 1) 


For  central  distributions  having  =  0;  1  <,j  <  r},  this  reduces  to  (a,  =  1/3;  1  <j  <  r}.  The 
resulting  transformation,  [(/„ ... ,  (/,]  —  [(6',/vp,)'/3 . (C/r/vp,)"3],  is  a  multivariate  exten¬ 

sion  of  the  Wilson- 1 Iilfertv  [32]  transformation  as  applied  to  each  variable.  Up  to  terms 
of  order  o(v~3),  expressions  for  the  means  [p,, ... ,  p,]  and  the  dispersion  parameters 
H  =  [£,,]  of  the  approximating  Gaussian  distribution  are  given  by 


Pj=l-2l9vpj,  1  <j<r 
ijj  =  2l9vPj;  \  <j<,r 


in  - 


2a2(y20  8<72((/20 


•-  *»»  ih. 


(3.12) 

(3.13) 

(3.14) 


The  corresponding  correlation  matrix  R  =  [p0]  has  elements 


o2(i/i)  4<r2(i/2/) 


Pij  „l/2„l/2  '  0  I  /2 

Pi  Pj  9vPt  Pj 


7 IT  [<t2((/20  -  (l/2)(p”+  )];  ‘*J-  (3-15) 
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Some  numerical  studies  are  reported  in  the  following  section  based  on  the  foregoing  results. 


4.  NUMERICAL  STUDIES 

Our  Gaussian  approximations  apply  to  various  joint  distributions  as  described  in 
Section  2.3.  Here  we  consider  implementing  these  approximations,  including  a  study  of 
their  accuracy  in  selected  cases. 

4.1  A  Monotone  Property.  Although  Gaussian  distributions  are  perhaps  best  known 
among  continuous  multivariate  distributions,  even  here  available  tables  are  limited  in  scope 
by  the  number  of  parameters  required.  The  following  facts  are  useful. 

Our  approximation  can  be  implemented  using  available  tables  and  a  result  of  Slepian 
[31]  which  asserts  that  if  PR( . )  and  Pr(  • )  are  Gaussian  measures  having  zero  means,  unit 
variances,  and  correlation  matrices  R  =  [p,,]  and  T  =  [y,J,  respectively,  such  that 
{y0<P,/,  1  <  i  <j  <  r),  then 

P /e(A'1  <  C|, ... ,  Xr  <  cr)  >  P p{A!  <  Cj, ... ,  Xr  <,  cr}  (4.1) 

for  any  scalars  (c,, ... ,  c,}.  In  the  special  case  that  [Z,, ... ,  Zr]  are  equicorrelated  standard 
Gaussian  variables  having  the  correlation  parameter  p,  the  cdf 

Fp(c . c)  =  />p(Z,  <  c . Zr<c )  (4.2) 

has  been  tabulated  in  [7]  for  various  values  of  r,  p  and  c.  If  we  now  take 
p*  =  min{  p,f,  1  <,  i<j  <  r],  Slepian's  inequality  assures  that  the  Gaussian  approximation 
for  [I7,,... ,  I7]  can  be  bounded  below  by  F,.(c„  — ,  cr).  In  the  central  case  for  which 
Q  =  0,  we  find  that  o2(ipi)  =  tr?,^,.  >  0  and  thus  p*  =  0  can  serve  as  a  lower  bound. 

4.2  Accuracy  of  the  Approximation.  We  compare  approximate  with  exact  probabili¬ 

ties  using  available  tables  or  algorithms  for  the  latter  in  selected  cases.  Two  cases  are  of 
interset,  namely,  (i)  the  bivariate  central  chi-squared  distributions  of  ( 20]  using  an  algorithm 
developed  in  (I5|,  and  (ii)  multivariate  central  chi-squared  distributions  studied  in  [22]  and 
tabulated  in  [21],  In  both  studies  the  approximating  Gaussian  probabilities  are  taken  di¬ 
rectly  from  [7],  To  avoid  interpolating  in  these  tables,  w'e  choose  values  of  r,  c  and  p  given 
there,  then  refer  to  (3.9)  to  express  the  inequalities  {Z,  <;  c;  l  <j  <,  r }  equivalently  as 
{U,£c*;  1  <,j<,r)  to  find  upper  limits  for  Ur~\.  Similarly  we  solve  expression 

(3.15)  with  p0  =  p  in  terms  of  p  to  find  the  parameters  appropriate  for  the  joint  distribution 
of  [£/,, ... ,  (/,].  Exact  probabilities  were  computed  in  the  fust  case  using  sufficient  terms 
of  a  series  given  in  [15],  and  in  the  second  case  by  interpolating  in  the  tables  of  [21|.  Details 
follow. 
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Let  v  =  l,r=  2,p=  2n,  Zn  =  Z22  =  L<  and  ?i2  =  r/„.  Then  the  joint  pdf  of 
{i\  =  trIFu(/2  =  trlF^}  is  given  in  [20|  as  a  series  bilinear  in  the  Laguerre  polynomials. 
An  algorithm  provided  in  [15]  was  programmed  and  used  to  compute  the  probabilities 
listed  in  Appendix  Table  B1  for  the  cases  {n=  1, 2,  3,  5,  15,  30  and  50}.  These  are  exact  to 
the  number  of  decimals  reported.  Four  conclusions  are  suggested,  other  factors  remaining 
constant,  (i)  For  any  n  and  t  the  error  of  approximation  tends  to  decrease  as  the  included 
probability  increases,  i.e.,  as  c  and  c *  increase.  In  connection  with  statistical  tests  having 
upper-tail  rejection  regions,  this  indicates  that  the  approximation  is  best  where  it  is  needed 
most,  (ii)  The  error  of  approximation  tends  to  decrease  as  n  increases.  This  is  indeed  the 
rationale  for  our  approximation,  which  is  here  seen  to  be  asymptotic  in  n  even  when 
v  =  1  in  our  earlier  developments,  (iii)  The  error  of  approximation  tends  to  increase  as  t 
and  thus  p  increase,  (iv)  The  Gaussian  approximations  appear  tenable  over  most  of  the 
parameter  values  studied,  especially  in  view  of  the  fact  that  accuracy  beyond  the  second 
decimal  is  seldom  required  in  practice. 

For  the  case  r  =  p  and  0  =  0,  the  distribution  of  [(/„  ... ,  f/r]  reduces  to  that  of  the 
diagonal  elements  of  a  central  Wishart  matrix  of  order  r  having  v  degrees  of  freedom  as  in 
[22|.  For  the  special  case  that  v  =  1  and  Z  =  [<r,J  is  a  correlation  matrix  with 
equicorrelation  parameter  6,  extensive  tables  of  this  distribution  are  given  in  [21J.  For 
(r=  3,4,6,  10}  and  various  values  for  the  parameters  p  and  c  of  the  approximating 
Gaussian  distributions,  actual  probabilities  and  their  Gaussian  approximations  are  given  in 
Table  B2  of  Appendix  B.  Exact  probabilities  were  obtained  from  (21)  using  a  standard 
two-dimensional  interpolation  procedure  ( cf.  |1],  p.  882)  to  interpolate  with  respect  to  c* 
and  S.  Several  conclusions  are  suggested  by  the  data  in  Table  B2  and  are  supported  by 
other  computations  not  reported  here.  Specifically,  the  error  of  approximation  tends  (i)  to 
increase  as  either  r  or  p  increases,  and  (ii)  to  diminish  towards  upper  tails  of  the  distrib¬ 
ution.  When  p  =  r  and  v  =  1,  we  accordingly  recommend  that  use  of  the  approximation 
be  restricted  to  upper  tails  only.  Despite  the  fact  that  our  developments  rest  heavily  on 
expressions  asymptotic  in  v,  our  approximations  fare  surprisingly  well  in  upper  tails  of  the 
distributions  studied  even  when  v  =  1 . 

APPENDIX  A 

To  supplement  developments  of  Section  3,  we  consider  first  the  joint  cumulants  and 
then  moments  of  [(/,, ... ,  (/,]. 

A.l  The  Joint  Cumulants.  Suppose  that  J?(W)  =  Z,  0).  The  joint  chf  of  IF 
with  symmetric  argument  7",  as  given  in  [2J,  yields  the  series  expansion  for  its  cgf  as  given 
in  [1 IJ  in  the  form 


II 


^CZ)“*S  is?-,s-,tr[crz)s-lr(z  +  Jg)]  (A.i) 

i=i 

where  £2  =  v-1©.  Standard  arguments  yield  the  joint  marginal  cgf  of  [t/,, ... ,  (/,]  on  setting 
other  arguments  to  zero  in  (A.  1),  giving  (3.1).  Letting  G{s)  =  (To^Y~lTo(^  +  *Q),  we  extract 
the  joint  cumulant  k  r  as  the  coefficient  of  i‘t[' ...  r^/s,!  ...  s,l  in  the  series 

*1/ (£)*=*£  rt-'s-'uGis).  {A. 2) 

S=  1 


Proceeding  recursively,  with  T0  and  Z  =  [Z0]  in  partitioned  form,  we  obtain 


(Tosr'=[z,  ■■  i , 


M.3) 


where  r'  =  s  —  1  and  the  expression  inside  brackets  is  the  typical  (u,v)  block  of  (T^L)’-'. 
Similarly  the  typical  block  of  G(s),  with  u  and  v  ranging  from  1  to  r,  is 


Quv(s)  ,  Z  ••• .  Z  (S4V  "b  sQi v) 

(2=1  1,  =1  * 

which,  on  summing  u  —  v  from  1  to  r,  yields 


(/1 .4) 


lrG(s)  =  tr[  dol)1-1  T0(L  +  r£J)] 

r  r 

=  V 


=  I  -  Z  -  \  tr[  &,t,  -  5vi,(5y,  + 

ii=i  i,=i 


M-5) 


Combining  (A. 5)  and  (A. 2)  yields  the  desired  cumulants  in  terms  of  the  parameters 
{Q,  =  tr(X,;  +  £!,,);  1  <j<r),  {co„.l  <j<r,s=  1,2,...},  and  Wi,  Wi  -  ViY  s  =  1.2,...}  as 
defined  in  Section  3.1.  These  are  summarized  in  Table  1  to  include  cumulants  of  order  s, 
with  1  ^  s  <  4.  Standard  relations  between  moments  and  cumulants,  as  listed  in  Table  Al, 
give  values  reported  in  Table  2. 


TABLE  Al.  Relationships  between  typical  central  and  noncentral  moments  and  cumulants 
of  order  s  =  r,  +  s2  +  r3,  1  <  s  <  4. 


s 

Relationship 

1 

Mioo  =  Kioo 

2  £  s  s  3 

=  K,l'2'3 

4 

=  Ktoo  + 

4 

M310  =  k310  +  2»CJ00K||0 

4 

M220  =  *220  "b  *200*020  b  2k}10 

4 

f*  211  =  *211  +  *200  *  011  +  2)CII0(C|01 
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A.2  Expansions  for  Moments  of  [Pj, ... ,  PJ.  Beginning  with  (3.4),  we  expand  terms 
appearing  there  in  their  binomial  series  to  obtain  expression  (3.5),  where  {a',  bf,  cj,}  are 
binomial  coefficients  given  by 

af  =  rx(rx  -  1) ...  ( rx  —  i  +  l)/t! 
b-  =sp{sp  -  i) ...  (r/?  -j  +  I)/ jl 
4  =  ty(ty  -  1) ...  (ty-k+  1  )/*! 

For  mixed  moments  of  order  two,  repeated  use  of  (3.5)  in  the  expression 
Mno =  Mno  —  Mioo  Moiu  gives 

S  I  (^-6) 

/= o  _/=o  '  7 


where  /((/,/)  =  [_E(e\e{]  —  E(e[)E(e{)'\.  Similarly,  from  the  identity  ^210  =  n'm  -  2/i'no  /^'00 
+  2(m!oo  )2  Foio  •  M200  Moio  and  (3-5),  we  determine  that 


m2,o(«./?)  =  .|  Z  (  V  )'+;  0Tl6?(c?-2atl)b}B{O,i,j) 
-2Z  I 


M-7) 


where  £(/t,  /,/)  =  £(ef)[£(e;ey)  -  E(e[)E(e{f).  Finally,  the  identity  Mm  =  Mm  -  mIoi  Moio  - 
Mill  Mi oo  -  Mho  Mooi  +  2mIoo  Moio  Mooi  after  lengthy  reduction  yields 


Mm(«,/?.y)  = 


oo  oo  oo 

m 


1=0  j= o  k=0 


(A.8) 


where 


Qi.j,  k)  =  [Efeje/e*)  -  E(e\e^)E(e{)  -  E(e[)E(e{e ■*) 
-  E(e\el)E{e^)  +  2£(ej)£(e2J')£(ej:)]. 


Expansions  (A. 6),  (A. 7)  and  (A. 8)  may  be  truncated  to  approximate  those  moments 
to  terms  of  specified  order  in  v.  It  is  required  to  compute  a  sufficient  number  of  terms  using 
expressions  for  E{e[e{e$)  from  Table  3  as  central  moments  for  [  i\ ,  U2,  6'J.  Observe  that 
{£(e,)  =  0;  1<;^3);  thus  expressions  for  ^,10(<*,  /?),  nm(x,f3)  and  nm(x,p,y)  simplify 
somewhat.  Specifically,  A(i,j)  =  0  =  Z?(0,  i,j)  if  either  i  =  0  or  j  =  0;  5(1,  i,f)  =  0  for  all  i 
and  j;  and  C(i,j,  k)  =  0  if  either  i—  0,  j  —  0  or  k  =  0.  Using  these  facts  and  properties  of 
the  binomial  coefficients,  we  summarize  in  Tables  4,  5  and  6  the  computations  needed  for 
Mno(®>  P)<  Mjio(“.  /?)>  and  Mm(«.  y)  vip  to  terms  of  order  o(v-3),  where  the  quantities  {£,,*} 
are  identified  in  Table  7. 
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APPENDIX  B 


TABLE  Bl.  Exact  probabilities  for  bivariate  x2  distributions  having  n  degrees  of  freedom, 
together  with  the  Gaussian  approximation*  (A)  having  parameters  c  and  p. 


p 

c 

1 

2 

3 

n 

5 

15 

30 

50 

A* 

0.1 

-2.0 

.0000 

.0001 

.0004 

.0005 

.0007 

.0008 

.0008 

.0009 

0.1 

-1.0 

.0197 

.0277 

.0290 

.0296 

.0303 

.0305 

.0307 

.0313 

0.1 

0.0 

.2736 

.2704 

.2688 

.2674 

.2663 

.2661 

.2660 

.2659 

0.1 

1.0 

.7154 

.7145 

.7144 

.7144 

.7145 

.7145 

.7145 

.7140 

0.1 

2.0 

.9542 

.9552 

.9555 

.9556 

.9556 

.9555 

.9555 

.9554 

0.1 

3.0 

.9977 

.9976 

.9976 

.9975 

.9974 

.9974 

.9973 

.9973 

0.2 

-2.0 

.0000 

.0002 

.0004 

.0007 

.0010 

.0011 

.0012 

.0014 

0.2 

-1.0 

.0215 

.0311 

.0331 

.0344 

.0360 

.0366 

.0369 

.0381 

0.2 

0.0 

.2897 

.2869 

.2853 

.2838 

.2825 

.2822 

.2821 

.2820 

0.2 

1.0 

.7276 

.7247 

.7238 

.7231 

.7223 

.7220 

.7218 

.7208 

0.2 

2.0 

.9558 

.9564 

.9565 

.9565 

.9563 

.9562 

.9561 

.9559 

0.2 

3.0 

.9977 

.9977 

.9976 

.9976 

.9974 

.9974 

.9974 

.9973 

0.4 

-2.0 

.0000 

.0002 

.0007 

.0012 

.0020 

.0023 

.0024 

.0029 

0.4 

-1.0 

.0257 

.0401 

.0440 

.0470 

.0501 

.0511 

.0517 

.0536 

0.4 

0.0 

.3236 

.3221 

.3202 

.3183 

.3163 

.3159 

.3157 

.3155 

0.4 

1.0 

.7489 

.7444 

.7426 

.7411 

.7392 

.7384 

.7380 

.7362 

0.4 

2.0 

.9590 

.9590 

.9589 

.9587 

.9583 

.9580 

.9579 

.9574 

0.4 

3.0 

.9978 

.9978 

.9977 

.9976 

.9975 

.9974 

.9974 

.9973 

0.6 

-2.0 

.0000 

.0003 

.0011 

.0023 

.0040 

.0045 

.0048 

.0055 

0.6 

-1.0 

.0321 

.0543 

.0606 

.0649 

.0687 

.0699 

.0705 

.0725 

0.6 

0.0 

.3623 

.3609 

.3585 

.3561 

.3535 

.3529 

.3527 

.3524 

0.6 

1.0 

.7690 

.7645 

.7626 

.7608 

.7585 

.7576 

.7571 

.7552 

0.6 

2.0 

.9624 

.9623 

.9620 

.9617 

.9611 

.9608 

.9606 

.9600 

0.6 

3.0 

.9980 

.9979 

.9978 

.9978 

.9976 

.9975 

.9975 

.9974 

0.7 

-2.0 

.0000 

.0004 

.0016 

.0034 

.0056 

.0063 

.0066 

.0074 

0.7 

-1.0 

.0372 

.0650 

.0723 

.0768 

.0804 

.0815 

.0820 

.0840 

0.7 

0.0 

.3844 

.3824 

.3798 

.3772 

.3746 

.3740 

.3737 

.3734 

0.7 

1.0 

.7794 

.7755 

.7737 

.7721 

.7699 

.7690 

.7685 

.7667 

0.7 

2.0 

.9643 

.9642 

.9640 

.9637 

.9630 

.9627 

.9625 

.9619 

0.7 

3.0 

.9981 

.9980 

.9979 

.9979 

.9977 

.9976 

.9976 

.9975 

♦Approximate  values  are  Fp(c,  c)  =  P„{ZX  <  c,  Z2  <  c)  from  |7], 
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TABLE  B2.  Exact  (E)  and  approximate*  (A)  probabilities  for  an  r-dimensional  x2  distrib¬ 
ution  having  v  =  1  degrees  of  freedom. 


0 

0.0 

1 

.0 

c 

2 

.0 

3.0 

r 

A 

E 

A 

IB* 

■B 

A 

r  = 

3 

0.1 

.1525 

.1489 

.6194 

.6106 

.9343 

.9343 

.9966 

.9960 

0.2 

.1741 

.1731 

.6451 

.6263 

.9383 

.9357 

.9967 

.9960 

0.4 

.2220 

.2232 

.6876 

.6597 

.9455 

.9393 

.9970 

.9961 

0.6 

.2810 

.2786 

.7249 

.6972 

A 

.9525 

.9459 

.9973 

.9963 

0.1 

.0868 

.0871 

.5413 

.5259 

r 

.9158 

.9140 

.9955 

.9946 

0.2 

.1077 

.1130 

.5794 

.5506 

.9226 

.9166 

.9957 

.9947 

0.4 

.1581 

.1692 

.6406 

.6007 

.9341 

.9233 

.9962 

.9948 

0.6 

.2267 

.2334 

.6922 

.6539 

.9447 

.9338 

.9967 

.9953 

0.1 

.0293 

.0331 

.4215 

.3972 

r  — 

O 

.8824 

.8756 

.9933 

.9920 

0.2 

.0434 

.0551 

.4791 

.4373 

.8953 

.8814 

.9938 

.9921 

0.4 

.0851 

.1100 

.5705 

.5144 

.9155 

.8958 

.9947 

.9925 

0.6 

.1569 

.1800 

.6450 

.5919 

.9325 

.9139 

.9958 

.9934 

r  = 

10 

0.1 

.0037 

.0066 

.2683 

.2404 

.8258 

.8065 

.9893 

.9867 

0.2 

.0079 

.0184 

.3469 

.2983 

.8515 

.S202 

.9903 

.9870 

0.4 

.0277 

.0605 

.4780 

.4072 

.8876 

.8507 

.9924 

.9881 

0.6 

.0822 

.1281 

.5853 

.5146 

.9154 

.8838 

.9943 

.9902 

♦Approximate  values  are  Fp(c, ... ,  c)  —  Pp(Z,  <  c, ...  ,Z,<,c)  from  j 7]. 
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